function [h,x,y] = lpcvd_mode2h_2D(z,L,method)
[m,n] = size(z);
assert((m-1)*2+1==n||(n-1)*2+1==m,'z must be square matrix');
mode = size(z,1)-1;
h = zeros(L,L);
dx = pi/(L);
x = 0:dx:(pi-dx);
y = 0:dx:(pi-dx);

if(nargin < 3)
    method = '';
end
% ----------------------------------------------------
% if (strcmp(method,'ifft') || strcmp(method,'fft'))
%     h = ifft2(z)*L^2/pi;
% else
%     for ix = 1:L
%         for iy = 1:L
%             for m = 0:mode-1
%                 for n = m:mode-1
%                     h(ix,iy) = h(ix,iy)+z(m+1,n+1)*exp(1i*2*(m*x(ix)+n*y(iy)));
%                 end
%             end
%             h(ix,iy) = h(ix,iy)/pi;
%         end
%     end
% end
% ----------------------------------------------------

% the structure of z
% z(m,n), 0<=m<=M,-M<=n<=M   z(-m+M+1,n+1);
n0 = mode+1;
if (strcmp(method,'ifft') || strcmp(method,'fft'))
    h = ifft2(z)*L^2/pi;
else
    for ix = 1:L
        for iy = 1:L
            for m = 0:mode
                for n = -mode:mode
                    if(m>0)
                        h(iy,ix) = h(iy,ix)+z(m+1,n+n0)*exp(1i*2*(m*x(ix)+n*y(iy)))+conj(z(m+1,n+n0))*exp(1i*2*(-m*x(ix)-n*y(iy)));
                    else
                        h(iy,ix) = h(iy,ix)+z(m+1,n+n0)*exp(1i*2*(m*x(ix)+n*y(iy)));
                    end
%                      h(iy,ix) = h(iy,ix)+z(m+1,n+n0)*exp(1i*2*(m*x(ix)+n*y(iy)));
                end
            end
            h(ix,iy) = h(ix,iy)/pi;
        end
    end
end
